Utiliser OpenStreetMap avec R - mise en pratique

Author

Louis Laurian, Ronan Ysebaert, Timothée Giraud

0. Packages

# 1. Packages qui interfacent des données OSM
# install.packages('maposm', repos = 'https://riatelab.r-universe.dev')
library(tidygeocoder) # géocodage
library(maptiles) # import de tuiles OSM (raster)
library(osmdata) # import de données OSM (vecteur)
library(maposm) # import de données OSM (couches géo)
library(osrm) # calcul d'itinéraires

# 2. Autres packages utilitaires
library(sf) # manipulation de données vectorielles
library(terra) # manipulation de données raster
library(mapsf) # cartographie thématique
library(maplegend) # légendes
library(spatstat) # analyse de semis de points
library(stplanr) # segmentiser des lignes
library(mapview) # cartographie interactive
library(lwgeom)

1. Géocodage

Le package tidygeocoder (Cambon et al. 2025) permet d’utiliser des services de géocodage en ligne. Par défaut, le package utilise Nominatim qui se base sur les données OSM.

pt <- data.frame(address = "28 Place de l'Hôtel de Ville, 60200 Compiègne")
pt <- geocode(.tbl = pt, address = "address", quiet = TRUE) |>
  st_as_sf(coords = c("long", "lat"), crs = 4326)
pt
Simple feature collection with 1 feature and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 2.826361 ymin: 49.41777 xmax: 2.826361 ymax: 49.41777
Geodetic CRS:  WGS 84
# A tibble: 1 × 2
  address                                                  geometry
* <chr>                                                 <POINT [°]>
1 28 Place de l'Hôtel de Ville, 60200 Compiègne (2.826361 49.41777)

2. Import de tuiles OSM

Le package maptiles (Giraud 2025) permet de télécharger des tuiles chez différents fournisseurs, dont OpenStreetMap, au format raster. 20 niveaux de zoom sont disponibles, selon la taille de l’objet géographique à cartographier. Un zoom de niveau 10 correspond à une échelle 1/500.000.

Plusieurs modèles de tuiles sont disponibles ici, certains nécessitent une inscription.

emprise <- st_buffer(pt, 3000) |>
  st_bbox()

tiles <- get_tiles(x = emprise,
                   project = FALSE,
                   crop = TRUE,
                   cachedir = "cache",
                   zoom = 14)

pt_reproj <- st_transform(pt, crs = "EPSG:3857")

mf_raster(tiles, expandBB = c(rep(-.2,4)))
mf_map(pt_reproj, pch = 24, cex = 1.3, col = "darkblue", lwd = 2, add = TRUE)
mf_annotation(pt_reproj, txt = "Mairie de Compiègne")
mf_title("Tuiles OSM - Défaut")

# Autre exemple de fournisseur
tiles <- get_tiles(x = emprise,
                    project = FALSE,
                    crop = TRUE,
                    provider = "OpenStreetMap.HOT",
                    zoom = 14)

mf_raster(tiles, expandBB = c(rep(-.2,4)))
mf_map(pt_reproj, pch = 24, cex = 1.3, col = "darkblue", lwd = 2, add = TRUE)
mf_annotation(pt_reproj, txt = "Mairie de Compiègne")
mf_title("Tuiles OSM - OpenStreetmap.HOT")

3. Import de couches géographiques OSM

Pour la cartographie, le package maposm (pas sur le CRAN). Il se base sur les fonctions du package osmdata (Padgham et al. 2023) (partie suivante) pour télécharger les données OSM, au format vectoriel.

La fonction principale, om_get(), permet de récupérer certaines couches géographiques utiles (routes, parcs, bâti, etc), autour d’un point, dans un rayon donné.

res <- om_get(x = c(st_coordinates(pt)[1],
                    st_coordinates(pt)[2]),
              r = 2000)
Getting urban areas: 0.84 sec elapsed
Getting buildings: 17.2 sec elapsed
Getting green areas: 1.58 sec elapsed
Getting roads: 0.47 sec elapsed
Getting streets: 1.33 sec elapsed
Getting railways: 0.48 sec elapsed
Getting water bodies: 2.78 sec elapsed
names(res)
[1] "zone"     "urban"    "building" "green"    "road"     "street"   "railway" 
[8] "water"   
mf_map(res$zone, col = "#f2efe9", border = NA, add = FALSE)
mf_map(res$urban, col = "#e0dfdf", border = "#e0dfdf", lwd = .5, add = TRUE)
mf_map(res$green, col = "#c8facc", border = "#c8facc", lwd = .5, add = TRUE)
mf_map(res$water, col = "#aad3df", border = "#aad3df", lwd = .5, add = TRUE)
mf_map(res$railway, col = "grey50", lty = 2, lwd = .2, add = TRUE)
mf_map(res$road, col = "white", border = "white", lwd = .5, add = TRUE)
mf_map(res$street, col = "white", border = "white", lwd = .5, add = TRUE)
mf_map(res$building, col = "#d9d0c9", border = "#c6bab1", lwd = .5, add = TRUE)
mf_map(pt_reproj, pch = 24, cex = 1.3, col = "darkblue", lwd = 2, add = TRUE)
mf_map(res$zone, col = NA, border = "#c6bab1", lwd = 4, add = TRUE)
mf_title("Autour de la mairie")
mf_arrow()
mf_scale(size = 500, scale_units = "m")
mf_credits("OpenStreetMap contributors, 2025", pos = "bottomleft")

4. Export de tags avec osmdata

Le package osmdata permet d’extraire des données OSM à partir de l’API Overpass. Il nécessite dans un premier temps d’établir une bounding box au sein de laquelle seront extraites les données, puis de l’exécuter (fonction add_osm_feature()). Le résultat est une liste, qu’il faut ensuite convertir en objet sf.

Les objets OSM prennent la forme de nodes (points), ways (lignes ou polygones), ou relations (relation entre plusieurs objets). Ces éléments sont décrits par un ensemble de tags, recensées sous forme de clé-valeurs (key-values). Ici, nous recherchons dans la clé amenity les points (nodes) correspondant aux restaurants, bars et cafés. En dehors des clés-valeurs principales, d’autres attributs peuvent être présents, comme en l’occurence le type de cuisine, ou bien les horaires.

# Définition d'une bounding box (emprise Compiègne)
q <- opq(bbox = emprise, osm_types = "node")

# Extraction des restaurants, bars et cafés
req <- add_osm_feature(opq = q,
                       key = 'amenity',
                       value = c("bar", "cafe", "restaurant"))

poi <- osmdata_sf(req)
poi <- poi$osm_points

# Sélectionner les variables
poi <- poi[, c("osm_id", "name", "amenity", "cuisine")]

# Intersection
poi <- st_intersection(poi, st_transform(res$zone, crs = 4326))

# Reprojection
poi_reproj  <- st_transform(poi , crs = "EPSG:3857")

La fonction om_map() affiche les couches géographiques récupérées.

# Affichage
om_map(res, title = "Autour de la mairie")
mf_map(poi_reproj, type = "typo", var = "amenity", cex = .7, pch = 15, 
       add = TRUE)
mf_map(pt_reproj, pch = 24, cex = 1.3, lwd = 2, col = "darkblue", add = TRUE)

Il est possible d’utiliser le package mapview (Appelhans et al. 2023) pour afficher les attributs de nos points, de manière interactive.

mapview(res$zone, alpha.regions = 0) + mapview(poi)
osmextract

Le package osmextract (Gilardi and Lovelace 2024) permet d’extraire des données directement depuis une base de données OSM, sans avoir à passer par une API. Le plus petit niveau géographique auquel il est possible de télécharger les données correspond en France aux anciennes régions. Pour travailler sur Compiègne, il nous faudrait d’abord télécharger l’ensemble des données OSM de la Picardie, disponibles ici.

5. Carroyage

Une fois les points d’intérêt récupérés, il est possible de les agréger dans une maille.

grid <- st_make_grid(res$zone, cellsize = 300, square = FALSE)
grid <- st_intersection(grid, res$zone)
grid <- st_sf(ID = 1:length(grid), geom = grid)

# compter
inter <- st_intersects(grid, poi_reproj, sparse = TRUE)
grid$n_poi <- lengths(inter)

om_map(res, title = "Autour de la mairie")
mf_map(grid, var = "n_poi", type = "choro", breaks = c(0,1,3,5,10,15),
       alpha = .6, border = NA, leg_title = "Nombre \nd'aménités",
       leg_val_rnd = 0, leg_pos = "topright", add = TRUE)
mf_title("Aménités gustatives - données carroyées")

6. Lissage (KDE)

Le package spatstat (Baddeley, Turner, and Rubak 2025) regroupe de nombreuses fonctions pour analyser des semis de points. Il utilise (génère ?) des objets de classe ppp (planar point pattern), dans des fenêtres d’observations, de classe owin (observation window).

La fonction density.ppp() calcule des valeurs à partir d’une estimation par noyau (KDE). Elle prend en entrée une portée (sigma), ainsi qu’une résolution au sein de laquelle sera calculée la densité, en pixels (eps).

Moraga (2023)
p <- as.ppp(X = st_coordinates(poi_reproj), W = as.owin(res$zone))
ds <- density.ppp(x = p, sigma = 100, eps = 10, positive = TRUE)

# Calcul densité de restaurants par hectare
r <- rast(ds) * 100 * 100
crs(r) <- st_crs(poi_reproj)$wkt
r
class       : SpatRaster 
dimensions  : 400, 400, 1  (nrow, ncol, nlyr)
resolution  : 10, 10  (x, y)
extent      : 312629.1, 316629.1, 6344049, 6348049  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / Pseudo-Mercator (EPSG:3857) 
source(s)   : memory
name        :         lyr.1 
min value   : 2.225074e-304 
max value   :  1.970817e+00 
mf_map(res$zone, col = "#f2efe9", border = NA)
mf_raster(r, leg_title = "N. restaurants / ha.", leg_pos = "left", 
          pal = c( "#75871B", "#BAB97D", "#FAF3CE", "#FDE16A", "#F9B40E",
                   "#E88C23", "#A25933"), add = TRUE)
mf_map(res$water, col = "#aad3df", border = "#aad3df", lwd = .5, add = TRUE)
mf_map(res$railway, col = "grey50", lty = 2, lwd = .2, add = TRUE)
mf_map(res$road, col = "white", border = "white", lwd = .5, add = TRUE)
mf_map(res$street, col = "white", border = "white", lwd = .5, add = TRUE)
mf_map(poi_reproj, pch = 20, cex = .2, col = "black", add = TRUE)
mf_map(res$zone, col = NA, border = "#c6bab1", lwd = 4, add = TRUE)
mf_title("La montagne alimentaire - Compiègne")

Il est possible de convertir cette couche en format vecteur. Par défaut, la fonction as.polygons() arrondit les valeurs des pixels à l’entier.

# par défaut
r_poly <- as.polygons(r) |>
  st_as_sf()

# tous les dixièmes
r_poly <- as.polygons(r, round = TRUE, digits = 1) |>
  st_as_sf()

mf_map(res$zone, col = "#f2efe9", border = NA)
mf_map(res$water, col = "#aad3df", border = "#aad3df", lwd = .5, add = TRUE)
mf_map(res$railway, col = "grey50", lty = 2, lwd = .2, add = TRUE)
mf_map(res$road, col = "white", border = "white", lwd = .5, add = TRUE)
mf_map(res$street, col = "white", border = "white", lwd = .5, add = TRUE)
mf_map(r_poly, var = "lyr.1", breaks = quantile(r_poly$lyr.1), border = NA,
       lwd = .2, type = "choro", leg_title = "Densité\n (n/ha)", alpha = 0.6,
       pal = "Mint", add = TRUE)
mf_map(poi_reproj, pch = 20, cex = .2, col = "black", add = TRUE)
mf_map(res$zone, col = NA, border = "#c6bab1", lwd = 4, add = TRUE)

Par aménité, discrétisation commune

# Cafés et bars ensemble
poi_reproj$amenity[poi_reproj$amenity %in% c("bar", "cafe")] <- "bar-cafe"
p <- list()
ds <- list()
r <- list()
sel <- levels(as.factor(poi_reproj$amenity))

for(i in 1:length(sel)){
  p[[i]] <- as.ppp(
    X = st_coordinates(poi_reproj[poi_reproj$amenity == sel[i] ,]),
    W = as.owin(res$zone))
  ds[[i]] <- density.ppp(x = p[[i]], sigma = 100, eps = 10, positive = TRUE)
  r[[i]] <- rast(ds[[i]]) * 100 * 100
  crs(r[[i]]) <- st_crs(poi_reproj)$wkt
  r[[i]] <- as.polygons(r[[i]], round = TRUE, digits = 1) |>
    st_as_sf()
}

# Récupérer l'amplitude
amp <- do.call(rbind, r)
amp <- sort(c(unique(amp$lyr.1)))
bks <- quantile(amp)
pal <- mf_get_pal(n = length(bks) -1, palette = "Reds", alpha = .6, rev = TRUE)

mf_map(res$zone, col = "#f2efe9", border = NA)
mf_map(res$water, col = "#aad3df", border = "#aad3df", lwd = .5, add = TRUE)
mf_map(res$railway, col = "grey50", lty = 2, lwd = .2, add = TRUE)
mf_map(res$road, col = "white", border = "white", lwd = .5, add = TRUE)
mf_map(res$street, col = "white", border = "white", lwd = .5, add = TRUE)
mf_map(r[[1]], var = "lyr.1", breaks = bks, type = "choro", border = NA,
       leg_title = "Densité\n (n/ha)", pal = pal,  add = TRUE)
mf_map(poi_reproj[poi_reproj$amenity == "bar-cafe" ,],
       pch = 20, cex = .2, col = "black", add = TRUE)
mf_map(res$zone, col = NA, border = "#c6bab1", lwd = 4, add = TRUE)
mf_title("Densité de bars et cafés")

mf_map(res$zone, col = "#f2efe9", border = NA)
mf_map(res$water, col = "#aad3df", border = "#aad3df", lwd = .5, add = TRUE)
mf_map(res$railway, col = "grey50", lty = 2, lwd = .2, add = TRUE)
mf_map(res$road, col = "white", border = "white", lwd = .5, add = TRUE)
mf_map(res$street, col = "white", border = "white", lwd = .5, add = TRUE)
mf_map(r[[2]], var = "lyr.1", type = "choro", border = NA, breaks = bks,
       leg_title = "Densité\n (n/ha)", pal = pal, add = TRUE)
mf_map(poi_reproj[poi_reproj$amenity == "restaurant" ,],
       pch = 20, cex = .2, col = "black", add = TRUE)
mf_map(res$zone, col = NA, border = "#c6bab1", lwd = 4, add = TRUE)
mf_title("Densité de restaurants")

Effet de bord

7. Temps d’accès

Le package osrm (Giraud 2024) interface l’engin de routage OSRM et propose différentes fonctions pour le calcul de temps de trajet et d’itinéraires. D’après la documentation de l’API OSRM, le nombre de requêtes est limité à 512 par connection à l’API, espacées d’au moins 5 secondes.

Pour des calculs plus volumineux (échelle nationale voire européenne), il est possible d’installer une instance OSRM sur son propre serveur (voir ici).

7.1. Depuis la mairie

Pour calculer les temps de trajet (en minutes) et les distances (en mètres) entre la mairie de Compiègne et l’ensemble des restaurants, bars et cafés situés dans un rayon de 2000m, nous utilisons la fonction osrmTable(), qui présente en sortie une liste composée de deux dataframes correspondant aux coordonnées des points d’origine et de destination, et d’un vecteur correspondant à la métrique calculée. Il est possible de calculer les distances et temps de trajet entre plusieurs points en une seule requête.

# Création de la matrice O/D
o <- data.frame(X = st_coordinates(pt)[, 1],
                Y = st_coordinates(pt)[, 2])
row.names(o) <- 1:nrow(o)

d <- data.frame(X = st_coordinates(poi)[, 1],
                Y = st_coordinates(poi)[, 2])

head(o)
         X        Y
1 2.826361 49.41777
head(d)
         X        Y
1 2.820433 49.41496
2 2.824740 49.41771
3 2.824230 49.41931
4 2.823624 49.42116
5 2.826944 49.41643
6 2.832764 49.41088
# 1 requête : 1 point d'origine, 92 points de destination
mat_dist <- osrmTable(src = o,
                      dst = d,
                      measure = "distance",
                      osrm.profile = 'foot')
Sys.sleep(5)

mat_dur <- osrmTable(src = o,
                     dst = d,
                     measure = "duration",
                     osrm.profile = 'foot')

dist <- t(mat_dist$distances)
dur <- t(mat_dur$durations)
poi$dist <- dist
poi$dur <- dur

apply(poi$dur, 2, summary)
                1
Min.     0.200000
1st Qu.  2.750000
Median   3.950000
Mean     6.087805
3rd Qu.  8.050000
Max.    22.700000
apply(poi$dist, 2, summary)
                1
Min.      16.0000
1st Qu.  204.0000
Median   295.0000
Mean     453.6098
3rd Qu.  596.0000
Max.    1693.0000

Il est possible d’atteindre 41 restaurants, cafés ou bars en moins de 4 min à pieds depuis la mairie de Compiègne, et les trois quarts se situent à moins de 596m à pieds.

7.2. Depuis la grille

nrow(grid) * nrow(d)
[1] 15498

Afin de construire des indicateurs d’accessibilité, il est nécessaire de calculer les temps de trajet vers tous les restaurants / bars / cafés de la ville depuis chaque point de grille.

Nous avons en tout 15498 couples origine / destination. Pour éviter de surcharger le serveur de démo, nous découpons nos requêtes en autant de destinations. Ainsi, chaque requête calculera le temps de trajet entre nos 189 points d’origine et un point de destination.

Les calculs ne seront pas présentés ici, et ont été effectués en local. Les résultats sont enregistrés dans le dossier data (fichier mat.rds).

osrmNearest

Pour que l’engin de routage calcule un itinéraire, les points d’origine et de destination doivent être localisés sur une route. Si le point d’entrée ne se situe par sur une route, il sera déplacé sur le point le plus proche situé sur une route, en distance euclidienne. Cela peut avoir un impact sur l’itinéraire final ainsi que les temps de trajet, plus ou moins important selon le type d’espace (urbain ou rural).

# Coordonnées des points d'origine : centroïdes des cellules de la grille
o <- st_centroid(grid) |>
  st_transform(crs = 4326)

o$X <- st_coordinates(o)[, 1]
o$Y <- st_coordinates(o)[, 2]
o <- st_set_geometry(o[, c("X", "Y")], NULL)

# Les points de destination restent les mêmes
mat <- list()

for(i in 1:nrow(d)){
  mat[[i]] <- osrmTable(src = o[, c("X", "Y")],
                        dst = d[i , c("X", "Y")],
                        measure = "duration",
                        osrm.profile = 'foot')
  Sys.sleep(5)
  mat[[i]]$ID_ORIG <- c(1:nrow(o))
  mat[[i]]$ID_DEST <- rep(i, nrow(o))
  mat[[i]] <- data.frame("ID_ORIG" = mat[[i]]$ID_ORIG,
                         "ID_DEST" = mat[[i]]$ID_DEST,
                         "duration" = mat[[i]]$durations)
  colnames(mat[[i]]) <- c("ID_ORIG", "ID_DEST", "duration")
}

mat <- do.call(rbind, mat)
head(mat)
saveRDS(mat, "data/mat.rds")

Nous pouvons maintenant calculer plusieurs indicateurs d’accessibilité : temps de trajet au restaurant le plus proche, au deuxième, au troisième (question du choix), nombre de restaurants à moins de 5 minutes, 10 minutes, 15 minutes, etc.

mat <- readRDS("data/mat.rds")
# Restaurant le plus proche pour chaque point de grille
poi_min <- aggregate(duration ~ ID_ORIG, mat, FUN = min, na.rm = TRUE)

# Récupérer l'identifiant de destination associé
poi_min <- do.call(rbind,
                     lapply(split(mat, mat$ID_ORIG),
                            function(x) x[which.min(x$duration), ]))

# Cartographier
grid <- merge(grid, poi_min,
              by.x = "ID", by.y = "ID_ORIG", all.x = TRUE)

om_map(res, title = "Temps d'accès au restaurant le plus proche")
mf_map(grid, var = "duration", type = "choro", breaks = seq(0,15,3),
       alpha = .6, border = NA, leg_title = "Temps\n(min.)",
       leg_val_rnd = 0, pal = "Greens", leg_pos = "topright", add = TRUE)

# Nombre de poi à moins de 5 minutes
mat5 <- mat[mat$duration < 5 ,]
mat5 <- aggregate(duration ~ ID_ORIG, mat5, FUN = length)

grid <- merge(grid[, "ID"], mat5,
              by.x = "ID", by.y = "ID_ORIG", all.x = TRUE)
grid$duration[is.na(grid$duration)] <- 0

om_map(res, title = "Restaurants à moins de 5 minutes")
mf_map(grid, var = "duration", type = "choro", breaks = c(0,1,5,10,30,40),
       alpha = .6, border = NA, leg_title = "Nb de\nrestaurants",
       leg_val_rnd = 0, pal = "Reds", leg_pos = "topright", add = TRUE)

8. Itinéraires

Pour aller plus loin, il est possible de récupérer les tracés des itinéraires, à l’aide de la fonction osrmRoute(). Attention, l’API route d’OSRM permet de calculer des itinéraires seulement entre deux points. Une requête correspond donc à un couple origine / destination.

routes <- list()
o <- data.frame(X = st_coordinates(pt)[, 1],
                Y = st_coordinates(pt)[, 2])
row.names(o) <- 1
for(i in 1:nrow(d)){
  routes[[i]] <- osrmRoute(src = o,
                           dst = d[i,],
                           overview = "full",
                           osrm.profile = 'foot')
   Sys.sleep(5)
}

routes <- do.call(rbind, routes)
saveRDS(routes, "data/routes.rds")

Le package stplanr (Lovelace, Ellison, and Morgan 2024) permet initialement de modéliser et de travailler sur des systèmes de transport. Nous l’utilisons ici à des fins cartographiques, pour découper des lignes en segments.

routes <- readRDS("data/routes.rds")
routes <- st_transform(routes, crs = "EPSG:3857")
routes$n <- 1
routes$distance <- routes$distance * 1000

# La fonction overline permet de compter le nombre d'occurence des tronçons
seg <- overline(routes, attrib = "n", fun = sum)

routes <- routes[order(routes$distance, decreasing = TRUE), ]
l_seg <- list()

# La fonction line_segment découpe un objet sf LINESTRING en segments réguliers
for(i in 1:nrow(routes)){
  l_seg[[i]] <- line_segment(
    routes[i, ],
    segment_length = max(routes$distance) / 20)
}

pal <- mf_get_pal(nrow(l_seg[[1]]), rev = TRUE, palette = "Blues 3")
om_map(res, title = "Routes", theme = "dark")
for(i in nrow(routes):1){
  mf_map(l_seg[[i]], col = pal, lwd = 1.6, add = TRUE)
}
leg(type = "cont", val = c(min(routes$distance), 1000, max(routes$distance)),
    pal = pal, pos = "left", title = "Distance (m)")

Le package osrm propose aussi des isochrones, à l’aide de la fonction osrmIsochrone().

o <- data.frame(X = st_coordinates(pt)[, 1],
                Y = st_coordinates(pt)[, 2])

isos <- osrmIsochrone(loc = o,
                      breaks = seq(0,30,5),
                      res = 60,
                      osrm.profile = "foot")

saveRDS(isos, "data/isos.rds")
isos <- readRDS("data/isos.rds")
isos <- st_transform(isos, crs = "EPSG:3857")
isos <- st_intersection(isos, res$zone)

om_map(res, title = "Isochrones autour de la mairie de Compiègne")
mf_map(isos, type = "typo", var = "isomax", alpha = .6, pal = "Geyser",
      border = "white", leg_pos = NA, add = TRUE)

Bonus : la tournée des pizzeria

Je me trouve à la mairie de Compiègne et souhaite goûter toutes les pizzas de la ville en un minimum de temps. Avec la fonction osrmTrip(), c’est possible !

Cette fonction répond au problème du voyageur de commerce, qui consiste à déterminer le plus court circuit passant par chaque point une seule fois.

# Coordonnées des pizzeria
pizz <- poi[!is.na(poi$cuisine) & poi$cuisine == "pizza" , "geometry"]

# Ajouter le point de départ : la mairie
pizz <- rbind(pt[1, "geometry"], pizz)
pizz$id <- 1:nrow(pizz)

trip <- osrmTrip(loc = pizz,
                 overview = "full",
                 osrm.profile = "foot")

saveRDS(trip, "data/trip.rds")
trip <- readRDS("data/trip.rds")
trips <- trip[[1]]$trip |>
  st_transform(crs = "EPSG:3857")

# Pour les labels : récupérer le début de la ligne
start <- lwgeom::st_startpoint(trips)
start <- do.call(rbind, start)
start <- data.frame(id = rownames(trips), X = start[,1], Y = start[,2])
start <- st_as_sf(start, coords = c("X", "Y"), crs = "EPSG:3857")

pizz <- poi[!is.na(poi$cuisine) & poi$cuisine == "pizza" ,]|>
  st_transform(crs = "EPSG:3857")

om_map(res, title = "La route de la pizza")
mf_map(res$zone, col = "#a83c0a33", alpha = .1, border = NA, add = TRUE)
mf_map(res$zone, col = NA, border = "#d7b578", lwd = 15, add = TRUE) 
mf_map(trips, col = "black", lwd = 2, add = TRUE)
mf_map(start[start$id != 1, ], pch = 20, cex = 2, col = "#496d0e", add = TRUE)
mf_map(pt_reproj, pch = 24, cex = 1.3, col = "darkblue", lwd = 2, add = TRUE)
mf_label(start, var = "id", cex = 1.1, overlap = TRUE, pos = 4, halo = TRUE)

References

Appelhans, Tim, Florian Detsch, Christoph Reudenbach, and Stefan Woellauer. 2023. Interactive Viewing of Spatial Data in r. https://CRAN.R-project.org/package=mapview.
Baddeley, Adrian, Rolf Turner, and Ege Rubak. 2025. Spatstat: Spatial Point Pattern Analysis, Model-Fitting, Simulation, Tests. https://CRAN.R-project.org/package=spatstat.
Cambon, Jesse, Diego Hernangómez, Christopher Belanger, Daniel Possenriede, and Otto Hansen. 2025. Tidygeocoder: An Intuitive Interface for Getting Data from Geocoding Services. https://CRAN.R-project.org/package=tidygeocoder.
Gilardi, Andrea, and Robin Lovelace. 2024. Osmextract: Download and Import Open Street Map Data Extracts. https://CRAN.R-project.org/package=osmextract.
Giraud, Timothée. 2024. Osrm: Interface Between r and the OpenStreetMap-Based Routing Service OSRM. https://CRAN.R-project.org/package=osrm.
———. 2025. Maptiles: Download and Display Map Tiles. https://CRAN.R-project.org/package=maptiles.
Lovelace, Robin, Richard Ellison, and Malcolm Morgan. 2024. Stplanr: Sustainable Transport Planning. https://CRAN.R-project.org/package=stplanr.
Moraga, Paula. 2023. “Spatial Statistics for Data Science: Theory and Practice with r.” https://www.paulamoraga.com/book-spatial/index.html.
Padgham, Mark, Bob Rudis, Robin Lovelace, Maëlle Salmon, and Joan Maspons. 2023. Osmdata: Import ’OpenStreetMap’ Data as Simple Features or Spatial Objects. https://CRAN.R-project.org/package=osmdata.